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ABSTRACT 


The analysis and incorporation into a multigrid scheme of 
several vectorizable algorithms are discussed. von Neumann 
analyses of vertical line, horizontal line, and alternating 
direction ZEBRA algorithms were performed; and the results were 
used to predict their multigrid damping rates. The algorithms 
were then successfully implemented in a transonic conservative 
full-potential computer program. The convergence acceleration 
effect of multiple grids is shown and the convergence rates of 
the vectorizable algorithms are compared to the convergence rates 
of standard successive line overrelaxation (SLOR) algorithms. 
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INTRODUCTION 


Rest. .* *. i 4 r the computation of transonic flows in recent 
years has been levoted to improving the accuracy, speed, and 
geometric capability of computational tools. Most three- 
dimensional transonic codes use the potential approximation to 
model the flow-field physics and the successive line overrelaxa- 
tion algorithm (SLOR) to solve the resulting system of equations. 
These three-dimensional transonic codes, however, still require 
large amounts of computer resources and are therefore expensive 
to use extensively. Thus, there is a need for improving the 
computational efficiency and reducing the cost of these codes. 

Two ways to increase the computational efficiency of a 
program are to improve its computation rate and the convergence 
rate of its algorithm. The computation rate of a program can be 
improved through efficient program coding and the use of faster 
computers (ref. 1). The efficiency of the coding of a program is 
programmer dependent; and it is not truly an area for research as 
much < s a topic for education. Improvements in computer speed 

are available bu - usually through specialized computer architect 

e 

ture such as the Control Data Corporation CYBER 203 and 
CYBER 205 computer systems. Vector architecture places certain 
restrictions on the algorithm in the program; and significant 
research has been devoted to the development of algorithms which 
can satisfy these restrictions. Algorithms which can make use of 
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the vector nature of the machines are referred to as vectoriz- 
able. Unfortunately, the workhorse algorithm of transonic 
potential codes, SLOR, is not completely vectorizable • Even if 
SLOR were completely vectorizable, it still has a very slow 
convergence rate for transonic flow problems* 

"■he ZEBRA algorithms of South, Keller, and Hafez (ref. 2) 
are a promising new class of vectorizable algorithms* Although 
vectorizable, the convergence rate for the ZEBRA algorithms is 
approximately the same as that for SLOR. The number of opera- 
tions required to obtain a solution is proportional to the square 
of the number of points (0(n 2 ) in the flow field, where n is 
the number of points in the flow field for potential flow cal- 
culations. Hence, for very fine grids, programs using the SLOR 
and ZEBRA algor' •' rs are too expensive for extensive use. 

The conet - of using multiple grids to accelerate the 
iterative solution of a set of finite difference equations was 
first proposed by Federenko (refs. 3 and 4) in 1961. The 
multiple grid concept was further analyzed by Bakhv^ov (ref. 5); 
but it was not until Brandt (ref, 6) ext> nded the technique and 
applied it to an elliptic set of equations that multigrid started 
to gain acceptance. 

Federenko showed that the solution for a set of equations 
could be obtained in 0(n) operations. Convergence in 0(n) 

operations using multigrid was proven for fairly general condi- 
tions by Nicolaides (ref. 7) and Hackbush (ref. 8). SLOR 

algorithms take 0(n 2 ) operations to obtain a solution, so it 
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is quite obvious that multiple grid (multigrid) techniques are 
quite attractive* 

In recent years , multigrid has gained wide acceptance as an 
acceleration technique* It was first applied to two-dimensional 
nonlifting transonic potential flow calculations in 1977 by South 
and Brandt (ref* 9) to accelerate an SLOR algorithm* This work 
was extended by Jameson (ref* 10) in 1979 to lifting airfoils but 
with an alternating direction implicit (ADI) type algorithm* 
Multigrid was applied to three-dimensional transonic potential 
flow calculations by Caughey (ref* 11) in 1983 using a penta- 
diagonal form of an SLOR-type algorithm* Although tremendous 
improvements in convergence rates were observed in each of these 
studies, the algorithms used with the multigrid limit the 
vectorizabili ty of the programs* 

In an attempt to take advantage of the fast convergence 
rate obtainable with multigrid and yet retain a high degree 
of vectorizability , multigrid is used to accelerate the 
vectorizable ZEBRA algorithms in the present study* Details of 
multigrid, the ZEBRA algorithms and their smoothing rates as 
compared to SLOR, and the results of the incorporation of ZEBRA 
algorithms in multigrid are given below* To reduce run costs, 
a two-dimensional problem is considered rather than a three- 
dimensional problem • 
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THE POTENTIAL MODEL FOR TRANSONIC FLOW 


Continuum Formulation 

Consider the flow of a perfect gas past an airfoil with a 
shape defined by the function Y(x). If irrotational flow is 
assumed, then the dimensional velocity, $, satisfies 

7 x $ » 0 


This condition is satisfied by the introduction of a reduced (or 
disturbance) potential function, 4>, such that 

V » v x + v $$ 

CO CD 

In conservation form, the two-dimensional potential equation is x 

(pu ) + (pv) * 0 ( 1 ) 

x y 

where p is the isentropic density: 


P 


00 


and 


a 2 « 1 /M 2 + (Y - 1) ( 1 - u 2 - v 2 )/ 2 

u = 1 + 4> x (2) 

v a <t> y 
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where p and a have been normalized by their free-stream values 
and u and v have been normalised by V a . Equation (1) is 
nonlinear and is strongly affected by the Mach number* This 
effect is more clearly seen if equation (1) is rewritten in 
nonconservation form: 

(a 2 - u 2 ) * xx + (a 2 - v 2 ) * yy - 2uv$ xy - 0 

Note that as u ♦ a (or M ♦ 1), the coefficient of $ xx 
approaches 0. If u > a (M > 1), the equation changes type and 
becomes hyperbolic instead of elliptic. 

In the present study, subcritical and supercritical flows 
over a symmetric, nonlifting 1 2-percent-thick parabolic arc 
airfoil are considered. The airfoil boundary condition is: 

3 Y v 

37 " u 

This boundary condition was applied at the airfoil chord line 
{ y = 0 ). Since only nonlifting cases are considered, symmetry 
boundary conditions ($ y = 0) are used fore and aft of the airfoil 
section. Zero-disturbance boundary conditions ($ = 0 ) are used 
at the upstream and downstream extents of the grid as well as the 
top boundary of the grid. (See fig. 1.) 
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Discrete Formulation 


Equation (1) may be solved at discrete points by rewriting 
it in a finite-difference formulation. The artificial density 
method is used in the present work to introduce dissipation. The 
discrete form of equation (1) with artificial density is t 

$(Fu) . 6(pv) 

» + — T = 0 

ox oy 


where 


M 0 - t lx P x , (3) 

P is the isentropic density, and 

t « max [ 0, 1 - (1/M) 2 ] (4) 

is the switching factor* 

In the present work, a modification to the switching factor, 
t f was used of the form: 

t * m ax [o # 1 - (M /M)^j (5) 

where M g is an input switch Mach number* By setting M g < 1 , 
the additional damping of the artificial density method can be 
introduced at subsonic points near the sonic region* This has 
been found to aid convergence (ref* 10)* 
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Boundary conditions as notad in figure 1 were used. Ghost 
points were used to allow central differences for the Neumann 
boundary conditions at the airfoil section and the symmetry line. 

A Cartesian grid with "interest fornula" grid stretching is 
used in the present study. In the y~dir*ction (see fig. 1). grid 
stretching begins at the airfoil where: 

Ay j+1 - P y * * yj (6) 

The parameter, F y is an input to the program and varies from 
1.0 for no grid stretching (uniform grid) to higher numbers such 
as 1.1 for 10-percent grid stretching. In the x-direction (see 
fig. 1), a uniform grid is used above the airfoil. Forward of 
the leading edge, 

&x i-l “ F x * Ax i (?) 
Aft of the trailing edge, 

Ax i+1 " F x * Ax i < 8 > 

The parameter, F x is an input to the program and is similar to 
Fy except it is used for the x-direction tretching. The local 
grid aspect ratio is defined as 

\ ■ Ax/Ay (9) 
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where Ax and Ay are the grid spacing in the x- and y- 
directions, respectively, as shown in figure 1. 

For the subcritical test cases run in the present study 
( = 0.1), a nonstretched uniform grid with X = 1.0 was used. 
In the supercritical cases (M* » 0.6), grid stretching was neces- 
sary to keep the outer boundaries far enough from the airfoil to 
eliminate their effect on the results. Therefore, 8 percent 
stretching was used in the x-direction (F x = 1.08) and 12 percent 
stretching was used in the y-direction (F y = 1.12). This gave a 
range of aspect ratios on the fine grid of 0.298 < X < 3.426. 

This same stretched grid was also used in the compressible sub- 
critical cases ( M a ■ 0.5). All cases were run using a 

fine grid with 65 points in the x-direction and 33 points in the 
y-direction unless otherwise noted. 

THE MULTIGRID METHOD 

The basis for multigrid is the use of successively coarser 
grids to calculate corrections to a fine grid solution. 
Excellent developments of the multigrid technique are given in 
references 4, 12, and 13. For completeness, brief developments 
of multigrid are given below for both linear and nonlinear 
operators . 
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Linear Equations 


Consider the problem 

L h U h = f h (10) 

where L h is a linear, finite-difference operator on a grid, G* 1 , 
with spacing h. The forcing function, f, is known and U h is 
the solution to the problem on the grid with spacing h. Xf we 
take u* 1 as an approximation to U* 1 with an error of 

V h » O h - u h , 

equation (10) can be written as 

L h (u h + V h ) = f h (11) 

Since L is a linear operator, this can be written as 

L h (u h ) + L h (V h ) « f h (12) 

Xf V 11 is a smooth function, it can be represented on a coarser 
grid, G^* 1 , with spacing 2h, twice the spacing between the points 
as the grid with spacing h. The grid G^* 1 is formed by 

deleting every other point in G* 1 . Therefore, G^€ G h . Points 
are eliminated from G^* 1 to form G^ and so forth to form 
G®* 1 , G 1 ®* 1 , etc. Bach subsequent grid is a subset of the previous 
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grid. In this context, a function is considered smooth if it 
contains no high-f requency components which will cause aliasing 
of the function when it is transferred to a coarser grid. 

It is possible to solve for an approximation to V* 1 on the 
coarser grid, G 2 * 1 , using equation (12) written for the coarser 
grid 

L 2h (l^ h v h ) » I 2h (f h - L h u h ) (13) 

IjJ is known as the restriction operator, and it simply transfers 

the values of a function from the fine grid to the coarse grid. 

Details of restriction operators are given later. For simplic- 

ity, define f 2h = I 2h (f h - L^u h ) as a forcing function on the 

coarse grid. Taking v 2 * 1 = I 2 **v h , equation (13) becomes 

n 

L 2h v 2h m f2h (14) 

Since equation (14) is for a coarser grid than equation (12), the 
numerical solution for V 2 * 1 is much cheaper to obtain because 
fewer points are involved. Or-* V 2h is obtained, it is used to 
correct the fine grid iterative solution, u h using 


(u h ) 


new 


(u n ) 


old 


+ I 


2h 


,2h 


(15) 
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The coarse grid to fine grid transfer operator, I*}. , in equa- 

2n 

tion (IS) is called the prolongation operator. This operator is 
discussed later. 

Since the form of equation (14) is identical to the form of 
equation (10), it is obvious that a grid with spacing 4h can be 
used to find corrections to the "solution" of the problem on the 
grid with spacing 2h. Successively coarser grids may be used 
until a grid is reached which is so coarse that the "solution" is 
obtained very, very quickly and cheaply. (The limit is a qrid so 
coarse that only a single unknown remains and is obtainable by 
direct solution.) The correction from the coarsest grid is then 
used to correct the correction on the next finer grid; and this 
is continued through successively finer levels until the finest 
grid is reached and the approximate solution is updated. 

The usefulness of corrections obtained on a coarser arid is 
dependent on the smoothness of the fine grid error passed to the 
coarse grid. Hence, it is absolutely necessary that the high- 
frequency components of the error on the fine grid are reduced, 
if not completely eliminated. It is the responsibility of the 
smoother (usually a relaxation algorithm) to damp the high- 
frequency components of the error. The low-frequency components 
of the error are unimportant for all but the coarsest grid since 
the grids which are coarser than the current grid will see these 
low-frequency components as high frequencies and they will be 
damped there. If the high-frequency components are not damped, 
then the restriction operator will pass aliased information to 



the coarser 


grid and the entire multigrid scheme will cease to 


converge* (An example of this is given in the Results 

section*) The choice of a smoother is critical to the use of 

multigrid* A subsequent section is devoted to a discussion of 
smoothers • 

The cycle of work performed starting on the finest grid# 
then visiting the coarser grids, and then returning to the finest 
grid is called one multigrid cycle* The cycles are repeated 

until sufficient convergence if obtained on the finest grid* 

In the present study, a fixed cycle known as the V-cycle was 
used* In the V-cycle, a prescribed number of iterations are 
performed on increasingly coarser grids, starting on the fine 
grid and proceeding to the coarsest grid, and then proceeding on 
increasingly finer grids back to the finest grid* 

Nonlinear Equations 

The previous development of the multigrid scheme was for 
linear operators* Since the potential equation used to model 
transonic flow is nonlinear, it is necessary to use the Full 

Approximation Storage (FAS) multigrid scheme (ref* 9) which is 
applicable to both linear and nonlinear problems* A development 
of the FAS scheme is given below in equations (16)-(20)* 

If the L* 1 operator is taken to be nonlinear, the step 

taken between equations (11) and (12) in the previous development 
is no longer valid* Instead, L h u h is subtracted from both sides 
of equation (11) to givet 
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For the coarse grid, this becomes: 



2hh 

n 





L h U h ) 


(17) 


If the second term on the left-hand side is moved to the right- 
hand side, equation (17) can be written as: 


L 2h (u 2h ) = f 2h 


(18) 


where 



Once values of are obtained, the fine grid iterative 

solution is updated using the following equation: 


(u ) 


new 


<u h ) 


old 


+ I 


2h 


[u 


2h 



( 20 ) 


Two points should be noted. First, the prolongated term on the 
right-hand side is a correction for the fine grid. Second, the 
operator used on the coarse grid (eq. (18)) has the same form as 
the fine grid operator, the grid spacing (h and 2h ) being the 


only difference 



Restriction Operators 


As mentioned previously# the restriction operator# » 

transfers the values of some discrete function from a fine grid 
to a coarser grid. The simplest restriction operator is 

injection. Here# the value of the function at each point on the 
coarse grid is equated to the value of the function at the 
coincident point on the underlying fine grid. In equation form 
for a generic function# 'I*# injection is: 


<l' 2h (x # y > 
o o 


'V y o> 


for all <x 0 , y Q ) C <G h n 0 2h ) . 

For some problems, it is necessary to use more complex 
restriction operators to transfer more "global" information from 
the fine grid. Normally# this is to eliminate small oscillations 
in the function on the fine grid. In the present work# a nine* 
point weighted average was used. 


<J» 2h (x # y ) - 1/4 <J> h (x , y ) 
o o o o 


+ 1/8 [tl» h (x + h, y ) + (x - h, y ) 
o o o o 


+ <l» h (x o # y o + h) + i|» h (x o # y Q - h)] 

♦ 1/16 [i|> h (x q + h, y Q + h) + i|» h (x o + h, y^ - h) 


+ V (x - h, y + h) + (x - h# y - h)J 
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2h 

Note that i|> (x , y ) is now influenced by the values of the 

o o 

function at nine points on G* 1 rather than just (x Q , y Q ) . 

Prolongation Operators 

L 

The prolongation operator, I^, the coarse- to fine-grid 
transfer operator. In the present study, direct transfer of the 
values of the function is used for fine grid points coincident 
with coarse grid points (Type A in fig. 2). Linear interpolation 
of coarse grid values is used for those fine grid points not co- 
incident with coarse grid points (Types B, C, and D in fig. 2). 
It is obvious that four steps are necessary to complete the 
prolongation . 


SMOOTHING ALGORITHMS 

As mentioned in the previous section, the smoother is used 
in a multigrid scheme to eliminate the high-frequency components 
of the error in the solution. The proper choice of a smoother is 
critical to the success of multigrid. Fortunately, it is usually 
possible to predict the performance of a smoother before it is 
incorporated in a multigrid context. 

Algorithms 

Six algorithms were considered as multigrid smoothers in the 
present work - vertical, horizontal, and alternating-direction 
SLOR and ZBBRA I. The three SLOR algorithms have been studied by 
other researchers and were included in the present work for 
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comparison with the ZEBRA I algorithms which are the focus of the 
current study* A description of each of the six algorithms is 
given below. 

VLOR - Vertical line overrelaxation* - When updating the 
point (i.j), VLOR uses updated values at three adjacent points; 
(i-1,j), (i.j+1). and (i,j-l), if the solution proceeds in 

th-% increasing i direction. (See fig. 3.) The values at 

I ' i ) are updated when the previous i « constant vertical 
li _s updated. The values at (i,j+1) and (i»j-l) are up- 
dated at the same time as the (i.j) point. The implicit set of 
tridiagonal equations generated by this scheme is readily solved 
using the Thomas algorithm. Since the solution of the implicit 
lines is order dependent and the Thomas algorithm is recursive. 
VLOR is not fully vectorizable • 

HLOR - Horizontal line overrelaxation. - The HLOR scheme is 
basically the same as the VLOR scheme except that the implicit 
lines are in the horizontal direction. When updating the 
point (i.j), HLOR uses updated values at three adjacent points; 
(i»j-l)» (i+l,j), and (i-i.j)» if the solution proceeds in the 

increasing j direction. (See fig. 4.) The values at (i.j-1) 
were updated when the previous j = constant, horizontal line was 
updated. The values at (i+l.j) and (i-1.j) are updated at 
the same time as the (i.j) point. The implicit set of tridiag- 
onal equations generated by this scheme is readily solved by 
the Thomas algorithm and. ae with VLOR. HLOR is not fully 
vectorizable . 
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ADLOR - Alternating direction line overrelaxation . - ADLOR is 
just the alternating application of VLOR and HLOR sweeps. Either 
VLOR or HLOR is used to update all of the points in the field and 
then the other is used to further update all of the points in the 
field. The ADLOR algorithm tends to be less sensitive to grid 
stretching. 

V2EB1 - Vertical line ZEBRA I. - VZEB 1 is a two-color, 
implicit, vertical line scheme. For the first color, say 
i « even points, only the "new" values at the po . its above and 
below and (i,j-i)) the point to be updated ((i,j)) are 

used. (See fig. 5.) This gives a set of tridiagonal equations 
which can be solved using the Thomas algorithm. Once the points 
in the first color are updated, their new values are used to 
update the points in the second color, odd values of i in this 
example. Hence, for the second color, updated values of the 
points (i-1,j), (i + t,j), (i,j— 1 ), and (i,j+1) are used to 
update the values at the point (i,j). (See fig. 5. ) The values 
at the points above and below the point at (i,j) and the values 
at (i,j) are updated simultaneously. This results in another 
set of tridiagonal equations which is solved using the Thomas 


algorithm. 

Since the implicit 

lines 

of a given 

color 

are 

decoupled 

from one another. 

all the lines of a given 

color 

may 

be solved 

simultaneously. 

Recall 

that 

the solutions of 

the 

tridiagonal lines in VLOR 

were 

order 

dependent so 

that 

the 

lines had 

to be solved one 

at a 

time . 

Hence, they 

were 

not 


vectorisable. The tridiagonal lines of a given color in VZEB1 
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are independent and can be solved simultaneously* This means 
that vector instructions can be used* Specif ically, vector in- 
structions are used to calculate the tridiagonal coefficients for 
all the implicit lines of a given color at each value of j* 
After all the coefficients are calculated, the back -subs titutions 
are performed for all the implicit lines, one value of j at a 
time* Further information about the vectorization o' the ZEBRA 
algorithms is contained in references 1 and 2* 

HZEBl - Horizontal line ZEBRA I* - RZEB1 is identical to 
VZEBl except that the implicit lines are in the horizontal 
direction • 

AD2EB1 - Alternating direction line ZEBRA I* - ADZEBi is just 
the alternate application of VZEBl and HZEBl sweeps* One of the 
two is used to update all of the points in the field and then the 
other is used to further update the points in the field* 

Smoothing Analysis 

The von Neumann stability analysis was developed by John 
von Neumann at Los Alamos around 1944* The method was circulated 
privately for several years (ref* 14)* A short description of 
the approach was first published in 1947 by Crank and Nicolson 
(ref* 15) and the first thorough explanation was given by 
O'Brien, Hyman, and Kaplan (ref* 16)* An excellent textbook 
explanation of the method was presented by Roache (ref* 17) in 
1 972* 
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In the von Neumann analysis, the approximate solution to a 
model equation generated by the algorithm of interest is expanded 
using a finite Fourier series • The exact solution is then 
subtracted - leaving the error* The decay or amplification of 
each frequency of the error is considered to determine stability 
or instability (ref* 17). The decay of the error is expressed by 
the von Neumann damping or amplification factor, g, as a function 
of frequency* The maximum value of g over the high-f reauency 
range of the error frequencies is defined as the smoothing factor 
(ref* 9), u* 

u 5 max {g(9 , 9 )}, */2 < |®„.» 9 | « 11 (2D 
x y x y 

where $ x and 0 y are the phase angles associated with the x- 
and y-coordinate directions, respectively* The full range of 
each 0 in a general solution is - it < 0 < it. The range from 
-x to x may be divided into two groups, the high frequencies 
from -x to -x/2 and x/2 to x and the low frequencies 
from -x/2 to x/2* The high-frequency components of the errors 
must be eliminated on the fine grid because they cannot be 
resolved on a coarser grid* The low-frequency components of the 
error are reduced on a coarser grid where one~half of them become 
"high" -frequency errors due to the change in grid spacing. 
Coarser and coarser grids are visited to eliminate lower 
frequency components of the rror as one-half of the low- 
frequency range of a given grid becomes the high-frequency range 
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on the next coarser grid. Therefore# the smoother only needs to 
eliminate the error frequencies from */2 < |$| < » and this is 
the range examined for the smoothing factor. (One additional 
restriction is that the error not grow in the low-frequency 
range, or g < 1 for 0 < |^| < ir/2.) 

A von Neumann analysis was performed for each of the six 
different algorithms considered in this study. In these 
analyses, generalisations for varying grid aspect ratio and Mach 
number were taken into account. Grid stretching was not directly 
modeled in the von Neumann analyses, but one of its effects was 
modeled by the inclusion of grid aspect ratio, X. 

Prom the description of the conservative transonic potential 
flow problem given previously, it is clear that Mach number 
effects must be included in the von Neumann analyses of the 
smoothers. Recall that the density appears as a coefficient in 
the conservative transonic potential equation. Variations In 
density due to compressibility effects cause nonlinearities. The 
type of the equation also changes with Mach number, elliptic 
for M < l and hyperbolic for M > 1. 

The conservative full-potential equation with artificial 
density is difficult to analyze using the von Neumann analysis. 
Considerably simpler is the analysis of the nonconservative 
small-disturbance potential equation (1 - M ) <t> xx + $ yy « 0. 
Results of the von Neumann analyses using the simpler equation 
are expected to be consistent with the full-potential equation 
due to the validity of the substitution as shown in Appendix A. 
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Hence, the performance of the smoothers operating on the 
transonic jmall-disturbance equation was analysed, rather than 
the conservative full-pctential equation. Two separate analyses 
were performed for each of the algorithms t One for M < 1 with 
central differencing and one for M > 1 with upwind differenc- 
ing. Over- or underrelaxation was accounted for in each subsonic 
analysis with the inclusion of the relaxation parameter, w. 
Additional draping in the form of B$ xt , where B is a free 
parameter, was included in the supersonic analyses. 

The von Neumann method for VLOR is straightforward and is 
included here to demonstrate the analysis for a simpler case 
before the more complicated details of VZEB1 are given. The 
analyses for the ZEBRA algorithms are more complex since they are 
two-color schemes. The development for HX»OR roughly follows the 
analysis of VLOR and so is not included. The analysis of HZEB1 
is given in Appendix B. The six algorithms analyzed and their 
amplification factors are summarized in Table 1. 

A von Neumann Analysis of VLOR 

The following analysis is for the nonconservative small- 
disturbance potential equation using central differencing in 
subsonic regions and upwind differencing in supersonic flow 
regions. Since the finite dif.'erencing operator is different for 
subcritical (central differencing) and supercritical (upwind 
differencing) flow regions, it is necessary to analyze the two 
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flow conditions separately 


The sabcritical case is analysed 


first. 

Subsonic Flow H < 1 

The smalx -disturbance equation in operator form is: 

[(1 - M 2 ) 5 + A 2 6 I » 0 for H < 1. 

' xx yy J ij 


where 

M is the local Mach number 

S xx is an undivided, central difference operator in the 
direction 

6yy is an undivided, central difference operator in the 
direction 

A is the grid aspect ratio, Ax/Ay 


and 




is the potential at the point i,j. 


VI* OR can be written as 


(1 


- [ E x 


- .n + 1 . , ~n+1 , _ 

A ^ + (- 21) + E 


V 1 

x^ijj 


A 2 (e + - 21 + E*) » 0 

V y yJ 


( 22 ) 


x- 

y- 


(23) 
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where 


E is a displacement operator in the x (E x ) or y (E y ) 
direction in the positive (E + ) or negative (E“) direction. 
That is. 


E (J) = 4 

X v ij y i*1 , j 


and 


E 


y ij 


i » j ±1 


is a value of the potential at the n** 1 time level 
previous to the update, and 4> 1 is an updated value of the 
potential which satisfies equation (23). 


Fourier coefficients are substituted for the error, e, which 
also satisfies the homogeneous difference equation given by 
equation (23). Therefore, 


e 


n 

ij 


n vC? (8+0 ) 

g e x y 


(24) 


Substitution of the Fourier components (eq. (24)) into equa- 
tion (23) for yields: 
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g 


n+1 

= 2 

n 

g 


, rr? 9 

(-1 )[ (1 - M 2 ) ( e X 


- * T) ♦ a- ==i ca cos e - 2 ) J 


2 o )-1 


(1 - M Z )(e 


-/TT 9 


( 25 ) 


- 2/w) + A /w (2 cos - 2) 


where g is the amplification factor; 0^ and 0^ are the phase 
angles associated with the x- and y-coordinate directions, 
respectively; and u> is the overrelaxation factor such that 




n + 1 




0) 


[? n+1 - * n ) 


Supersonic Flow M > 1 

For a Mach number greater than or equal to 1*0, the small- 
disturbance potential equation is* 


r 

i 


(i 


M 2 ) $ + A 2 6 - 86 J * 

xx yy xt J ij 


0 for M > 1 


(26) 


where 


M, A, 6 , and <p . . 

yy ij 


are as defined below equation (22) 
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is an undivided upwind difference operator in the x- 
direction 

6 is a free parameter 


and 


6 ^ is an undivided second difference operator in the x- 
and time-directions which is added for stability. The 
difference expression is: 


6 . 
xt y i} 


U n+1 - * n ) . (* n+1 - a n ) 

^ 9 ij 9 i j ' 9 i-l,j J 


Using the operators defined with equation (23), equation (26) may 
be written as: 


(1 - M 2 ) ( (E‘) 2 - 2E‘ + i) + X 2 (E _ 

v x x * ij y 


21 + E + ) <J> 

y 


n + 1 
ij 


+ 


8 [[I - E") (J. 


n+1 

ij 



( 2 ?) 


Substitution of the Fourier coefficients and simplification 
yields : 
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g b B[l - (cos 0^ - sin e x > ]/{ (1 - m 2 )[(cos 9 x - / -1 ' sin 9 X ) 2 

- 2 (cos ® x - / -1' sin 9 x ) + l] + X 2 (2 cos 9^ - 2) 

+ 3 [l - (cos 9^ - / -1 sin 9^)]} (28) 

Prom equations (27) and (28) it is clear that VLOR is a marching 
scheme in purely supersonic flow and, if 6 » o, g * 0 for all 
frequencies* 

For a given method, Mach number (M), and values of 8, w, 
and X; values of g may be found as a function of 9 X and 
9 y where 9 X and 9 y vary from -180° to 180°. It is useful to 
look at contour plots of g for varying 9 X and 9 y as shown 
in figure 6. The frequency range for each of the 8's can be 


divided into two segments - 

the 

high-frequency 

range 

(-1 80° 

to 

-90° and 90° to 180°) and 

the 

low-frequency 

range 

(-90° 

to 


90°). Mote that the maximum values of g in figure 6 are in the 
low-frequency ranges of 9 X and 9 y (g max “ 0.97). 

With multigrid, only the high-frequency region is important 
and the low-frequency region may be ignored as reflected in 
figure 7. The maximum value of g in the high-frequency range 
is the smoothi.ng factor, u • Notice that the elimination of the 
low-frequency region significantly lowers the maximum value of 
g obtainable and so u is significantly less than 9 ma x ^ or 
this case (w a 0.44 versus 9 ma x “ 0.97). 
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A von Neumann Analysis of VZEB1 
Subsonic Flow M < 1 

Using the previously defined displacement operators, the 
VZEBl algorithm applied to equation (22) yields: 


2 , - n ~n+1/2 + n . 

(1 - M ) IE <(> - 21$ + E ()> 1 

x ij ij x ij 


2 r - +, ~n+1/2 

+■ X [E - 21 + e J $ 
y y i j 


r "n+1/2 n-1/2 - , n 

- B i $ - $ - E f <J> 

* ij ij x ij 


n- 1 . , 

♦u 11 


( 29) 


where : 


<t> 


n 

ij 8 


a new value of the potential at the "other" color 
at the n tfl time level. The two colors are 
considered to be one-half time step apart in the 
analysis . 


$ 


n-1/2. 

ij 


an old value o:: the potential on the color being 
updated . 


1 ) 


a new value of the potential on the color being 
updated (overrelaxed). 


$ 


n + 1/2 

i j / 


a value of the potential on the color being 
updated which makes the residual zero. 
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The three time levels in equation (29) are related by 


❖ 


n + 1/2 

ij 


1/ID 0 


n + 1/2 

ij 


u-1 ,n-1/2 
— *ij 


( 30 ) 


where u> is the relaxation parameter* Substitution of the 
Fourier coefficients for the potential and simplification gives a 
quadratic equation for /g 1 . 


2 2 2 
[(1 - M )(-2)(1/u>) + X (2 cos 9 - 2)(1/<o)J (/g) 

V 


[ ( 1 - M ) ( 2 cos 6 ) j /g’ 

x 


r 2 . CO-1 , 2 , <0-1 . . 

+ [(1 - M ) ( -2 ) ( ) + X (2 cos 0 - 2) ( )J - 0 {31) 

to y <0 


This may be solved using the quadratic equation to find 

the two roots* The root which represents g in the physical 

problem cannot be determined. The purpose of the von Neumann 

analysis is to predict the multigrid smoothing rate for 

VZEB1 . Since the value of u is governed by the maximum value 

of g , it is prudent to pick the root of g which has the 

maximum absolute value* 

Supersonic Plow M > 1 

In operator form* VZEB1 applied to equation (26) can be 
written as t 
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(1 


2 

- M ) 


f f E "> 2 

i >. * > 

X 


n-1/2 

ij 


2 (E ) 

x 


n n+1/2 

+ U> 

ij 


2 r - +, n+1/2 

+ X [E + 21 + E J <t> 

y y i j 


. n + 1/2 n - \ / 2 . - w n n-1 

• 8 l*i, - ♦„ - E X J - *11 11 


( 32 ) 


There are some interesting things to note here* First, an 
"old" value is used at i-2 to maintain the independence of the 
vertical implicit lines. Second, note that the terms represent- 
ing 6 are split between the two time steps. This equation 

xt 

may be rewritten as a cubic equation in /g*. 

a + b ( /g) 2 + c /"g* + d » 0 (33) 


where 

a “ (1 “ M 2 ) + X 2 (2 cos 0^ • 2) * 3 

b = [(1 - M 2 )(-2) + Bj (cos £ - /^7 sin 9 ) 

X X 

e ■ (1 - M ) (cos 6 - /-7 sin 0 ) + B 

d * - ft ( cos 0 /IT sin 0 ) 

X X 
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The maximum magnitude of the three roots 


is taken as the 


magnitude of g. 

In figure 8(a), a contour plot of g for VZEB1 at the 

conditions of M w ■ 0 • 1 , 3 ■ 0, X ■ 1.0, and u ■ 1,0 is given. 

These parameter values are the same as figures 6 and 7 for 
VLOR. Note that there are three areas of g * 1.0 in this 
contour plot. There is one at O x , ©y) equal to (0°, 0°) as with 
VLOR and there are two at (-180°, 0°) and ( + 180°, 0°). The 

latter two areas of g * 1 are not in the low-frequency regions 
and so are not discarded when u is examined, figure 8(b). For- 
tunately, since the values of 0 equal to 0° and ±180° are not 

obtainable on a discrete grid, only a limited effect of this 
region should be felt in a multigrid application. However, 
since u for VZEB1 is 0.97 and u for VLOR is 0.44 for the 
given conditions, it is anticipated that the VLOR may give a 
better convergence rate as a multigrid smoother. 

It is useful to consider only u rather than the full 

range of values of g. In figure 9, u for VZEB1 is plotted 

versus Mach number for three values of aspect ratio. X • 1.0 
is the uniform grid aspect ratio used for the entire grid for 
the incompressible test case mentioned below in the Results 
section. X « 3.426 and 0.0298 are the upper and lower limits of 
aspect ratio on the fine grid of the stretched grid in the 
compressible subsonic and supersonic test cases also discussed 
in the Results section. It is clear from figure 9 that the 
performance of VZEB1 is strongly affected by aspect ratio over 
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the entire Mach number range presented* The predicted perform- 
ance is poor for X < 1 but better for X > i, as expected for 
VZEBl . 

A von Neumann analysis for HZEB1 is given in the Appendix B* 
Contour plots of g for HZEB1 are quite similar to those of 
VZEB1 except, that they are rotated 90°* This is reflected in 
figure 10 where u versus M is plotted for various values of 
aspect ratio* Note that HZEB1 has good performance for small 
aspect ratio and poor performance at high aspect ratio, the 
opposite of VZEB1* 

The von Neumann analysis of ADZEB1 is just a combination 
of the analyses of VZEB1 and HZEB1* At each frequency, the 
damping rate for ADZEB1 is the square roo: of -h « product of the 
individual damping rates of VZEBl and HZEB1* This reduces the 
sensitivity of ADZEB1 to aspect ratio effects sfnce frequencies 
that either VZEBl or HZEB1 has problems with are compensated for 
by good performance by the opposite smoother (HZEB1 or VZEBl, 
respectively; i*e*, for small-aspect-ratio qrid cell , HZEB1 does 
well and for large-aspect-ratio grid cells, VZEBl does well* The 
improved performance of ADZEB1 for varying aspect ratio is 
reflected in figure 1 • 
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RESULTS 


A computer code was written to test the performance of the 
ZEBRA algorithms as multigrid smoothers* The transonic full- 
potential equation was modeled using central differencing with 
the artif* 1 density method used to provide upwind bias in 
supersonic egions. Small-disturbance boundary conditions were 
used to model the presence of an airfoil. 

The flow over a symmetric* 1 2-percent-thick parabolic arc 
airfoil was calculated. Free-stream Mach numbers of 0.1, 0.5, 

and 0.8 were considered. Stretched grids (8 percent in the x- 
direction and 12 percent in the y-direction) were used at all 
three Mach numbers. A uniform grid (no stretching, Ax * Ay) was 
also used at M s ■ 0.1. Only nonlifting cases were considered. 

The performance of VZEBl was examined at the conditions of 
Mg, s o.i on a uniform unstretched grid. Simple injection of the 
potential and residual was used in a 5-grid multigrid scheme. 
The solution did not converge. To understand v'.y this occurred, 
surface plots of the potential, residual, and forcing function 
for the various grids were made. 

A surface plot of the potential on the fine grid (grid 1 ) 
just before injection to grid 2 is shown in figure 12. The 
coordinate system for the grid is as shown in the figure and will 
apply to all other surface plots. Note that the positive 
direction for the function shown is down. Note also that the 
function is reasonably smooth. 
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A surface plot of the residual on grid 1 is shown in 
figure 13. Note the jagged nature of the residual function. If 
the x-y plane is rotated up and toward the viewer., so that the 
line of sight is directly down the y-axis (see fig. 14), it can 
be seen that every other "x - constant" line has a residual of 
sero. This is to be expected with the decoupled-line, two-color 
nature of the ZEBRA I schemes. 

The high-frequency oscillations in the residual on the first 
grid caused aliasing of the forcing function on the second 
grid. (The forcing function, f, is the difference of the fine 
and coarse grid operators; see eq. (19).) With no residual 
weighting (see Restriction Operator section above), the forcing 
function on the second grid had the shape shown in figure 15. 
With residual weighting, the forcing function on the second grid 
had the shape shown in figure 16. 

The effect of residual weighting on the convergence history 
is shown in figure 17 where the logarithm of the ratio of the 
fine grid residual to the initial fine grid residual is plotted 
versus work. One unit of work is defined as the number of 
operations required to perform one fine grid (G* 1 ) iteration. 
One iteration on grid G^* 1 requires 0.25 work units for a two- 
dimensional problem, etc. 

From the above discussion, it is obvious that residual 
weighting must be used with VZEB1 to eliminate aliasing of the 
coarse grid forcing function and allow convergence. (Further 
research showed that, for t'.ese same reasons, residual weighting 
is necessary for HZEB1 also.) 
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Weighted averaging of the potential function was also 
examined. Zt was found that this was not necessary since the 
potential is relatively smooth on the fine grid (fig. 12) and the 
injected values of the potential on the coarse grid are not 
aliased without the weighting. The "shape" of the coarse grid 
potential is shown in the surface plot of figure 18. This same 
shape is obtained independent of potential weighting. Injection 
was used for the potential in all subsequent results discussed. 

with residual weighting, VZBB1 was greatly accelerated by 

the use of multigrid. In the most benign case considered in 
this study (M n e 0.1, uniform grid), the acceleration of VZEB1 
by multigrid is illustrated in figure 19 using 2, 3, 4, and 

5 grids. All runs were conducted with an overrelaxation factor 
of 1.0. Note that 4- and 5-grid multigrid give identical 

convergence and are only slightly faster than 3-grid multigrid 
for these conditions. At a Mach number of 0.1, the flow is 
completely subsonic and effectively incompressible. The 

potential equation governing the flow is elliptic; and, since 

the aspect ratio is unity, it effectively reduces to Laplaces' 

equation. Thus, the excellent speedup obtained using multigrid 

was expected even though it was better than the von Neumann 

analysis suggests. Evidently, the frequencies giving the 
smoothing factor (u) as plotted in figure 9 were not obtained. 

Figure 8(b) shows that the maximum values of g for VZEB1 at 
M = 0.1 and X ■ 1.0 are concentrated at the two points of 

high frequency in 9 X and the low frequency in 9 y . If 
these frequencies do not appear in the error of the solution. 
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then they do not need to be damped and the effective u is 
reduced. This evidently occurred at * 0.1 and X » 1.0 

for VZEB1. Prom this example, it is clear that the results from 
a von Neumann analysis must be carefully examined to determine 
validity. Predictions of u can be biased by frequencies not 
found in the coded problem. 

The performance of VZEB1 multigrid is dependent on the 
grid aspect ratio. This dependence is reflected in figure 20. 
Starting with the most benign conditions, » 0.1 and an un- 

stretched grid, the performance of VZEB1 multigrid is excel- 
lent. If the Mach number is held the same but grid stretching 
is included (8 percent in the x-direction and 12 percent in 
the y-direction) the performance is significantly degraded. 
(Compare M® « 0.1 curves.) Note that not only are the absolute 
levels of the residual higher but also the asymptotic rate of 
convergence is not nearly as good with the grid stretching. 
Recall that the von Neumann analysis did predict a sensitivity to 
aspect ratio. 

If the Mach number is increased slightly to M m = 0.5 so 
that compressibility effects are in the solution, the performance 
is slightly worse than the incompressible case on the same 
stretched grid. Most of the difference is in the absolute levels 
of the residual and not in the asymptotic convergence rate. A 
further increase in Mach number to Mq» * 0.8 so that the flow 
field becomes transonic further degrades the performance of VZEB1 
multigrid but again mostly in absolute levels of residual and not 
asymptotic rate of convergence. (Compare stretched grid curves 
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in f ig • 20.) 


From this discussion it is clear that the 


performance of VZEB1 multigrid is affected more strongly by qrid 
aspect ratio than by compressibility (nonlinearity) or super- 
critical flow regions (equation type changes). This behavior was 
supported by the von Neumann analysis as summarized in figure 9. 

The acceleration afforded VZEB1 by various numbers of grids 
in the multigrid scheme at « 0.8 is shown in figure 21* The 

speedup at » 0.8 is not as dramatic as at » 0.1. Note 

that 4- and 5-grid multigrid give nearly identical performance 
and are only slightly better than the 3-grid multigrid* The 
asymptotic rate for 5 grid is only slightly better than for 
1 grid, but the absolute level of residual is significantly less* 

Since many others have used VLOR as a multigrid smoother, it 
is worthwhile to compare the performance of VZEB1 and VLOR in 
multigrid* At - 0*8 on a stretched grid, their performance 

is compared in figure 22* Note that their rate of convergence is 
nearly identical, however the VZEB1 can be more fully vectorized 
and so should give better overall performance (less computer time 
to obtain a converged solution)* Similar comparison plots for 
HZEB1 multigrid are shown in figures 23-26* The results follow 
those for VZEB1* It should be noted that for the test problem in 
the present work, HZEB1 gives better performance than VZEB1* 
This is not a general result, but is obtained here only because 
of the specific grid used and the difference in percentage of 
total points included in the update implicit lines of the two 
algorithms • 
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In the von Neumann analysis of ADZEB 1 , it was found that the 
alternating direction scheme should give better performance on 
varying aspect ratio grids. (See fig* 11*) The prediction is 
confirmed in figure 27 where the effects of M w and aspect ratio 
on ADZEB 1 multigrid are examined. Note that the incompressible 
(unstretched and stretched grid) and the compressible cases all 
give roughly the same convergence rates. With the effects of 
grid stretching reduced, the effect of supersonic regions is seen 
to adversely affect the convergence rate. The performance of 
ADZEB 1 in the transonic flow calculation is quite good and is 
compared to ADLOR in figure 28. Unfortunately, because of the 
alternating direction nature of ADZEB1, vectorization is not as 
straightforward as with VZEB1 and HZEB1. 
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CONCLUSIONS 


Three vectorizable multigrid smoothers were examined; 
vertical, horizontal, and alternating direction ZEBRA. The 
smoothing performance of each was predicted using von Neumann 
stability analyses. The effects of Mach number, grid aspect 
ratio and damping factor were included* As expected, the 
analyses predicted that VZEBI performs best on grids with X > 1 
and HZEB1 performs best on grids with X < 1 . The analysis of 
ADZEB1 predicted that it is less sensitive to grid aspect ratio 
than VZEBi or HZEB1 • All three ZEBRA methods were predicted to 
have poorer performance at supercritical Mach numbers than at 
subcritical Mach numbers. 

The actual performance of each of the ZEBRA algorithms was 
then assessed by incorporation into a two-dimensional full- 
potential code. The performance of the vectorizable algorithms 
was compared to three well-known SLOR algorithms. In each 
comparison (VZEBI to VLOR, HZEB1 to HLOR, and ADZEB1 to ADLOR ) , 
the ZEBRA scheme had a convergence rate comparable with the 
respective SLOR scheme. These comparisons indicated a reasonable 
level of success of the multigrid acceleration of each ZEBRA 
scheme. Using very powerful smoothers, it is possible to have 
faster multigrid convergence rates than were obtained in the 
present study* However, the purpose of the present work was to 
study vectorizable ZEBRA multigrid algorithms, so slightly less 
than optimal convergence rates were acceptable as a compromise to 
obtain vectorizable code* 
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It is unfortunate that the von Neumann analyses did not 


predict quantitatively the performance of the ZEBRA schemes. In 
a discrete problem, only a finite number of frequencies are 
present in the error. These frequencies are dictated, in part, 
by the grid. A von Neumann analysis models a continuum problem 
and so considers all frequencies of the error. The value of u 
obtained from the analysis can be biased if it is predicted based 
on frequencies which are not present in the discrete problem, 
particularly when 0 V or equals or is close to 0 or in, 

a y 

Finally, it should be noted that only a nonlifting airfoil 
was considered in the present study. Convergence rates of cal- 
culations for lifting airfoils using the current methods would 
probably not be as fast. 
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APPENDIX A 


SIMPLIFICATION OF CONSERVATIVE FULL-POTENTIAL 
EQUATION WITH ARTIFICIAL DENSITY 


Consider the one-dimensional form of the full-potential 
equation using the artificial density method: 


(pu) « 0 (A1 ) 

x 


where 


P = p - tAxp 

X 


( A2 ) 


P 



a 2 } 1/(Y 


1 ) 


(A3) 


and 


® 2 - 1/ M 2 + (Y - D [1 - u 2 )/2 


(A4) 


(Note that equations (A1)-(A4) are the same as equations (1)-(4) 
in the main paper with v = 0* ) The finite difference form of 
equation (A1 ) is : 


P i+1/2 U i+1/2 " P i-1/2 U i-1/2 
Ax 


R ± (4« ) 


( A5 ) 


where 
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p i+1/2 ° P i+1/2 “ T i (p i+l/2 " p i-1/2 ) 

, , ( A6 ) 

" (1 " T i P i+1/2 + i P i-1/2 

and t is the switching factor as described in equation (5). 

A 

The residual, R^(4>), can be expanded about the exact solution , 
such that R^ ($ ) ® 0, 

* i + 1 3 R^ 

r" (<>) = R i («t>) + l UJJ - * k ) 


(A7) 


i + 1 3R" „ 

- . r in 

0 * *.L ^ Sk 


where e£ is the error in the solution at the n* 1 * 1 time level* 
Substituting equation (A6) into equation (A5) and taking the 
partial derivative of Ax 2 R^, we obtain: 


6(Ax 2 R i ) = (1 - T) uAx«P i+1/2 + TAxu5 P i .i /2 + pAx6u i+ i/2 


- (1 - x) uAx6p 1 _ 1/2 - TAx « 6 P i . 3/2 " pAx6u i-1/2 


" (p i+1/2 ‘ P i-1/2 ) Ax6 (T i u i+1 /2 ) 


* (p i-1/2 " p i-3/2 J Ax6 (T i-1 U i-1 /2 ) 


(A8) 


41 



Since ( p i+1/2 “ p i-1/2* Ax ar *d * p i-1/2 ~ p i-3/2^ Ax are 

O(Ax^), the last two terms in equation (A8) may be dropped 
because they are higher order terms than the first six terms. 
The values of the derivatives in equation (A8) are: 


dp 

, 1 + 1/2 _ „2 
uAx — tr L — = pM 

’’i + l 


3p 


U A x 


ill LI B p„- 


3 


9fl l-1/2 


- pM 


»Ax ^-V 2 - PH 2 
H i-1 


dp. 


uAx 


77 


i -3/2 


- pM 


i-1 


A 3p i-3/2 
uAx 

*♦ i-2 


PM 


{ A9 ) 


PAX . p 




i + i 


du 


pAx 


i + 1/2 

TT 


- p 


.. 3 °i -1 / 2 

pAx — ■ p 


du 


pAx 




i-1/2 


- P 


i-1 


Substituting equations (A9) into equations (A7) and (Ae) and 
simplifying gives: 


K 2 « n 

Ax R^ 


- t 1 - (1 - T ) M 2 ] ej +1 1(2 - 3 x) M 2 - 2 ) ej 


t ( 1 - 3 t) m 2 - 1 ] ej -1 - T M 2 e ”_ 2 


(A10) 
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For M < 1, equation (4) gives T ■ 0. Substitution of this into 
equation (A10) gives: 


Ax 2 r" « p (1 - M 2 ) [ej +1 - 2eJ + «i +1 J (All) 

For M > 1, equation (4) gives t ■ 1 - 1/ M ^« Substitution of 

this into equation (A10) gives: 

Ax 2 rJ = P (1 - M 2 ) [ej - 2e”_ 1 + e"_ 2 J (A12) 

From equations (All) and (A12) we see that the conservative full- 

potential equation with artificial density may be simplified to 
the nonconservative small-disturbance potential equation with 
upwind differencing in the supercritical regions for the purposes 
of analysis.* 


*The development showing this simplification is based on work 
originally done by Jerry C. South, Jr., and is contained in his 
personal notes. 
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APPENDIX B 


ANALYSIS OF HZEB1 


Subsonic Flow - M < 1 

The small! 'disturbance equation in operator form is 


[ ( 1 - M 2 ) 6 + X 2 6 I ♦ . . - 0 

1 xx yy J ij 


( B 1 ) 


where each term is defined after equation (22) in the main body 
of this text. U?ing the displacement operators defined in the 
text, the HZEB1 algorithm applied to equation (B1) yields: 


(1 - M 2 ) [b“ - 21 + E + j + X 2 (-21) 

1 x x J ij ij 


+ X 2 ( E* + E + ) » 0 

• Y y J ij 


( B2 ) 


The time levels of the potential are defined with equation (29) 
in the text. Substitution of the Fourier coefficients for the 
potential and simplification give a quadratic equation for 
This can be solved for the two roots of /g* which are squared 
to get two roots of g. Which of the roots represents g in 
the physical problem can.'* be determined. The purpose of the 
von Neumann analysis is to predict the multigrid smoothing rate, 
v, for HZEB1. Since u is governed by the maximum value of g, 
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it is prudent to pick the root of g from the quadratic equation 
which has the maximum absolute value* The two roots are: 


g 


l'g) 


2 




2 


( B 3 ) 


where : 


1 [(1 - M 2 ) (2 cos e - 2) - 2X 2 ] 

U) 1 X J 


b ® 21 cos 9 


and 


c a (U) - 1 ) a 


Supersonic Flow - M > 1 

HZEB1 applied to the transonic small-disturbance equation 
for supersonic flow can be written as: 

!" - " 2 > *t-J A *ij - <»«> 

where 


n,n+i/2 
ij ' 


«*> 


«T 

o 

XV 



+ \ 4 


(♦ 


n + 1 /2 
ij +1 



+ 


n+1 /2 
ij-1 


) 
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and the time levels of 1> are defined above* The inclusion of 

"new" values at i - 2 would result in a pentadiagonal set of 
equations. To eliminate this, a $ xt type term is added. For 
stability, additional <fr xt is also added. The following 4 xt 
terms are used: 


- <1 - 


M 2 ) 


( A 


ij 


' & »i-2, 3 > 


+ S (A* 


ij 


* A »i-2, 3 > 


Adding these terms to the left-hand side of equation (B4) gives: 


r -2 (1 - M 2 ) + S ! (A*. - A£ . 

ij 


) + X A$ 


ij 


i, . » 2 , 




(B5) 


Substitution of the Fourier coefficients for the potential 
followed by algebraic manipulation of the equation results in a 
quadratic equation for /g*. 


a 



+ b » g + c * 


0 


(B6) 


where : 


— / - 10 

a * X 2 + [-2 (1 - M 2 ) + Sj (1 - e X ) 


b = - 2 X cos 9 
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and 


A 2 + 


(1 - 


M 2 ) 


+ (1 


M 2 ) 


- 2/-1 


9 

x 


Although the above represents a complete von Neumann 
analysis, HZEB1 was studied in greater detail in the present 
work* Values of g can easily be found from equation (B6), 

but it is not immediately obvious that the algorithm is stable, 
g < 1 for all (9 X , 9 y ). It is possible to show this 

analytically. Equation (B4) can be rewritten in the form: 

2 

z + nz + Y * 0 (B7) 


where 

z = /g’ 
n = b/a 

and 

Y = c/a 

Since n and y are complex, the conditions that must be 
satisfied for z (and so g) to be less than or equal to 1 are: 

Y | < 1 

and (B8) 

n - yn 1 < i - | y | 2 
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where n is the complex conjugate of n. The inequalities in 
equation 1 B 6 ) lead to the restriction: 


3 > 2 - 2M 4 


/ 


( 1 - M) ( 1 - b - 


M 2 , 


for M > 1 


( B9 ) 


A numerical experiment was performed to determine the sharpness 
of the restriction in equation (B7). A slightly supersonic free- 
stream case (M,,, » 1 . 01 ) was run with a very thin airfoil ( 0 . 10 - 
percent thick parabolic arc) to produce a nearly uniform super- 
sonic flow. Various percentages of the "minimum" 6 prescribed 
in equation (B7) were used as the coefficient of the additional 
9 xt term. It was found that the 3 required by equation (B9) 
was not exact. As little as 80 percent of the 6 required 

by equation (B9) was found to provide stability. Although the 
von Neumann analysis does not exactly predict the performance 
of HZEB1 , it does provide a qualitative indication of the 
performance . 
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Figure 1 . - Coordinate system and schematic of stretched grid* 
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Figure 2.- Prolongation operator- 
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Figure 3.- Schematic of VLOR algorithm* 
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Figure 4.- Schematic of HLOR algorithm. 
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Figure 5.- Schematic of VZEB1 algorithm* 
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von Neumonn Domping Role Contours 
Equation: (1 - M 2 ) * 0 

Method: VLOR 


Porometers: 
u = 1.000 

A = 1.000 

M = .1000 

fi « 0 . 
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Figure 6*- von Neumann damping rate contours for VLOR. 
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Multigrid Smoothing Factor Contours 
Equation: (1 - M 2 ) 4> xx + ^ = 0 

Method: VLOR 

Parameters: 

u = 1.000 

A = 1.000 

M = .1000 

fi * 0 . 


Figure 7.- Multigrid smoothing factor contours for VLOR. 
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Parameters: 
u = 1.000 
A = 1 .000 

M « .1000 

0 = 0 . 


(b) Multigrid smoothing factor contours. 


Figure 8.- von Neumann damping rate and multigrid smoothing 
factor contours for VZEB1 . 
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Figure 9 


\ = 3.426 0 

= 1.000 = 

= .2980E-01 O 



M 

Equation: (1 - M 2 ) o xx + <t> yy -• 0 

Method: VZEB1 

u = 1.000 

0 - 1 .000 


Multigrid smoothing factor versus M and X for 
VZEB1 . 
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Figure 10.- Multigrid smoothing factor versus M and 
HZEB1. 
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Figure 11.- Multigrid smoothing factor versus M and X for 
ADZEB1 . 
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Figure 14.- End view of residual on finest grid usin • ' !EB 1 

( Mgo « 0.1). 
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MAXIMUM |f | = 0.4323 



Figure 16.- Forcing function on second finest grid using VZEBl 
and residual weighting (M* * 0.1)* 




Work 


Figure 17.- Effect of residual weighting on VZEB1 (M w = 0*1)* 
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Figure 18 . « R 
( 



Figure .i.- Effect of multigrid on VZEB1 («„, « 0.1 )♦ 
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Figure 20. - Effect of M w and grid stretching on VZEB1. 
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figure 23.- Effect of multigrid on H2EB1 (M^ = 0.1). 
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Figure 26.- Comparison of HZEB1 and RLOR in muXtigrid (M w «* 0.8, 
5 grids). 


C - J. 
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Figure 27.- Effect of M*, and grid stretching on AOZBBi. 
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Log R/R 



Figure 28.- Comparison of ADZEB1 and ADLOR in multigrid (M w = 0.8. 
5 grids). 
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